Computational drug design of novel COVID-19 inhibitor

Background In 2003, the first case of severe acute respiratory syndrome coronavirus (SARS-CoV) was recorded. Coronaviruses (CoVs) have caused a major outbreak of human fatal pneumonia. Currently, there is no specific drug or treatment for diseases caused by SARS CoV 2. Computational approach that adopts dynamic models is widely accepted as indispensable tool in drug design but yet to be exploited in covid-19 in Zaria, Nigeria. In this study, steps were taken to advance on the successful achievements in the field of covid-19 drug, with the aid of in silico drug design technique, to create novel inhibitor drug candidates with better activity. In this study, one thousand human immunodeficiency virus (HIV1) antiviral chemical compounds from www.bindingBD.org were docked on the SARS CoV 2 main protease protein data bank identification number 6XBH (PDB ID: 6XBH) and the molecular docking score were ranked in order to identify the compounds with the highest inhibitory effects, and easy selection for future studies. Results The docking studies showed some interesting results. Inhibitors with Index numbers 331, 741, and 819 had the highest binding affinity. Similarly, inhibitors with Index number 441, 847, and 46 had the lowest hydrogen bond energy. Inhibitor with index number 331 was reported with the lowest value (− 48.38kCal/mol). Five new compounds were designed from the selected six (6) compounds with the best binding score giving a total of thirty (30) novel compounds. The low binding energy of inhibitor with index no. 847b is unique, as most of the interaction energies are of H-bond type with amino acids (Thr26, Gly143, Ser144, Cys145, Glu166, Gln189, Hie164, Met49, Thr26, Thr25, Thr190, Asn142, Met165) resulting in an overall negative value (−16.31 kCal/mol) making it the best of all the newly designed inhibitors. Conclusions The novel inhibitor is 2-(2-(5-amino-2-((((3-aminobenzyl)oxy)carbonyl)amino)-5-oxopentanamido)-4-(2-(tert-butyl)-4-oxo-4-(pentan-3-ylamino) butanamido)-3-hydroxybutyl) benzoic acid. The improvement it has over the parent inhibitor is from the primary amine group attached to meta position of first benzene ring and the carboxyl group attached to the ortho position of the second benzene ring. The molecular dynamics studies also show that the novel inhibitor remains stable after the study. This result makes it a better drug candidate against SARS CoV 2 main protease when compared with the co-crystallized inhibitor or any of the 1000 docked inhibitors. Supplementary Information The online version contains supplementary material available at 10.1186/s42269-022-00892-z.


Background
Coronaviruses (CoVs) have caused a major outbreak of human fatal pneumonia since the beginning of the twenty-first century. COVID-19 is an acute respiratory disease caused by the RNA virus SARS-CoV-2. In severe cases, the infection can cause pneumonia, severe acute respiratory syndrome, kidney failure, and even death (Du et al. 2020). There is currently no specific medicine or treatment for diseases caused by SARS-CoV-2 (Huang et al. 2020).
SARS-CoV-2 virus targets cells through the viral structural spike (S) protein that binds to the angiotensinconverting enzyme 2 (ACE2) receptor. Once inside the cell, viral polyproteins are synthesized that encode for the replicase-transcriptase complex. Structural proteins are synthesized leading to completion of assembly and release of viral particles (De Groot et al. 2013;Lau et al. 2005;Reusken et al. 2013).
Hosseini and Amanlou (2020) conducted a virtual screening procedure employing docking of 1615 FDA approved drugs to identify new potential small molecule inhibitors for protease protein of COVID-19 and their result indicates that among all FDA-approved drugs, simeprevir which is used for Hepatitis C virus (HCV) NS3/4A protease inhibitor, revealed strong interaction with protease binding pocket and placed well into the pocket even better than the lopinavir-ritonavir (Abd El-Aal et al. 2022;Al-Hossainy et al. 2021;El Azab et al. 2021). Since this compound is FDA-approved and has successfully passed various testing steps, they suggested that this drug could be a potential drug for treating the COVID-19 (Hosseini & Amanlou 2020). Motiwale et al. (2020) and friends applied molecular docking approach in conjugation with molecular dynamics (MD) simulations to find out potential inhibitors against Mpro of SARS-CoV-2 from previously reported SARS-3CL protease inhibitors. They used a total of 61, previously known inhibitors, where 4- (Sacco et al. 2020) benzoic acid, and 4-(4-methoxyphenyl)-6-oxo-2-[(2-phenylethyl)sulfanyl]-1,6-dihydropyrimidine-5-carbonitrile were reported to have minimum and maximum binding energy, respectively (Motiwale et al. 2020).
To achieve a fast and reliable drug in this current crisis, we initiated a virtual screening procedure (Gagic et al. 2020), employing docking of 1000 HIV1 protease inhibitors compounds from www. bindi ngBD. org over binding pocket of SARS CoV 2 main protease using 1 pdb file PDB: (6XBH) downloaded from Research Collaboratory for Structural Bioinformatics (RCSB) which represent main protease of SARS CoV 2 to identify potent inhibitors against the virus and to design a novel drug whose molecular dynamics studies will be done to ascertain the effectiveness in inhibiting the SARS-CoV-2 Mpro (Sacco et al. 2020).
The molecular docking result would provide first-hand knowledge about the interactions between the ligands and the target receptors since most of these ligands work by profoundly inhibiting the specificity and efficiency of protein (target) action (Adeniji et al. 2020;Stockwell 2000). All these will undoubtedly offer important structural insight into the design of novel COVID-19 drugs (Abd El-Aal et al. 2022;Al-Hossainy et al. 2021;El Azab et al. 2021). Besides, these methods will help to: promote savings in the cost of drug design and development, reduce the requirement for lengthy and expensive animal tests and, promote green chemistry to increase efficiency and eliminate chemical waste (DiMasi et al. 2003). Overall, the purpose of the work is to carry out a computeraided drug design (CADD) of SARS CoV 2 main protease inhibitor.

Ligand preparation
The 2D structure of each inhibitor was drawn using the ChemDraw v16.0 Windows 10 (32 bit and 64 bit), Copyright 1998-2016 PerkinElmer Informatics Inc and Page 3 of 29 Arthur et al. Bulletin of the National Research Centre (2022) 46:210 presented in table (Arthur et al. 2020). The structures were introduced into wavefunction 14 graphic user interphase (GUI) after which the 2D structures were converted into 3D structures by selecting the view dialog box present on Spartan 14 GUI. From the build option on Spartan 14, the structures were clean by checking to minimize using molecular mechanic force field (MM+) option to remove all strain from the molecular structure. In addition, this will ensure a well-defined conformer relationship among compounds of the study (Viswanadhan, Ghose, Revankar, & Robins, 1989). From the setup calculation option on Spartan 14, the calculation was set to equilibrium geometry at the ground state using a semiempirical PM6 (Parameterization Method 6) (Bikadi and Hazai 2009).
There are five different interaction potentials that contribute to the overall free binding energy established between the receptor pocket and the docked ligand (Gallicchio et al. 2010). These potentials include van der Waals potential for a hydrogen atom probe, van der Waals potential for a heavy-atom probe (generic carbon of 1.7A radius), hydrophobic energy terms, optimized electrostatic energy term, and lone-pair-based potential, which reflects directional preferences in hydrogen bonding. These energy terms are based on the all-atom vacuum force field with added functions to account for solvation free energy, desolvation energy and entropic contribution. It was shown that after each random step, full local minimization greatly improves the efficiency of the procedure (Abd El-Aal et al. 2022;Al-Hossainy et al. 2021;Ibrahim et al. 2020;Mohamed et al. 2022;Zwawi et al. 2021). The ICM program relies on global optimization of the entire flexible ligand in the receptor field and combines large-scale random moves of several types with gradient local minimization and a search history mechanism (Arthur et al. 2018).
The protein structure was surrounded with a 15 Å layer of TIP4P BOX water molecules. The electrostatic charge was neutralized by adding counter ions using the LeaP program of AMBER ver.11. After minimization, heating and equilibration, the production MD phase was carried out at 300 K for 1 ns with a time step of 1 ps (picoseconds) using the constant volume and temperature (NVT) ensemble and the Particle Mesh Ewald algorithm for the calculation of electrostatic interactions (Haspel, Zheng, Aleman, Zanuy, & Nussinov, 2017). The initial velocity of atoms was generated at 100 K in heating phase with a Maxwellian distribution and maintained. The pressure was kept at 1 bar by Berendsen weak coupling approach during equilibration (Berendsen et al. 1984).

Results
The numerical results of this study are presented in the tables presented below. This has become necessary because of the need to correlate some of the data. Other results, such as plots and pictorial representation of the interactions between the ligands and their receptor binding sites are presented as figures. Table S2 shows the molecular docking result of the reference inhibitor and 1000 HIV 1 antiviral compounds on SARS-CoV-2 main protease receptor (PDB ID: 6XBH). The following parameters are shown: number of rotatable torsions, hydrogen bond energy, hydrophobic energy in exposing a surface to water, the van der Waals interaction energy (sum of gc and gh van der Waals), internal conformational energy of the ligand, the desolvation of exposed H-bond donors and acceptors, the solvation electrostatics energy change upon binding and mean force Page 4 of 29 Arthur et al. Bulletin of the National Research Centre (2022) 46:210 score. According to the molecular docking results, it was found that the binding energy of co-crystallized ligand, (R)-3-(((2R,5S)-5-(((S)-(benzyloxy)(hydroxy)methyl) amino)-1-hydroxy-4-oxo-6-phenylhexan-2-yl) amino)-1,3-dihydro-2H-pyrrol-2-one was − 23.56kCal/mol while the binding energy of all the 1000 HIV 1 antiviral inhibitors lies between − 4.73 and − 48.38 kCal/mol. Figure Table 1. The ICM score for the best possible interaction pose was given as − 23.56 kCal/mol. The molecular docking score was ranked in order to identify the compounds with the highest inhibitory effects, and easy selection for future studies. The docking studies showed that inhibitors with Index numbers 331, 741, and 819 had the highest binding energies of all the compounds that were docked on SARS CoV 2 main protease (Additional file 1: Table S2). Similarly, inhibitors with Index number 441, 847, and 46 had the highest hydrogen bond energy. Inhibitor with index number 331 was reported with the highest value (− 48.38 kCal/ mol) and inhibitor with index number 46 having the least value (− 15.69 kCal/mol) for all the best six (6) selected. The high correlation of H-bonds with the number of flexible bonds (nflex) reflects on the high binding energies of the compounds, with the exception of inhibitor with index number 331 which has six (6) flexible bonds. The low binding energy of inhibitor with index number 331 is unique, its amplified hydrogen bond energy was as a result of an inductive effect created by the presence of three Π-sulfur interactions observed with Cys145, Met165 and Cys145, an amide pi stacked interactions with Thr24 and Thr25, π-lone pair interaction with Thr24, π-pi stacked interaction with Hie41 and π-alkyl interaction involving Met49 and Cys145. Another noticeable point is the bond length of two of the conventional H-bonds involving Thr26 and Hie164 with bond length 1.68 and 1.94A, respectively, thus impacting positively on its activity.
Based on binding energy ranking, inhibitors with Index numbers 331, 741, and 819 were selected for design. Similarly, inhibitors with Index numbers 441, 847, and 46 were selected based on hydrogen bond energy for design. Five new inhibitors labeled a-e were designed for each of the selected inhibitors above. All the compounds in the dataset docked were found to inhibit the receptor by completely occupying the active sites in the target receptor. Most of the inhibitors were tangled in both hydrophobic and hydrogen bonding interactions with the receptor. Here, it was found that strong inhibitor binding is reflected by the frequency of hydrogen bonds. In addition, the molecular docking studies carried out show that all the compounds were found to inhibit the receptors by completely occupying the active sites in the target receptor. The mechanism for this reaction is the same in all cases, which includes the intercalation of the inhibitors between the covalently bonded SARS CoV 2 main protease complex. Additional file 1: Table S3 shows the structures and IUPAC name of designed novel inhibitors while Table 8  Other noticeable interactions with the receptor include π-alkyl interaction mediated through Cys145. The inhibitor benzyl (5-amino-1-((4-(2-(tert-butyl)-4-oxo-4-(pentan-3-ylamino) butanamido)-3-hydroxy-1-phenylbutan-2-yl) amino) -1,5-dioxopentan-2-yl) carbamate (Index number 847) from which it was designed has binding score energy of − 39.89 kCal/mol and H-bond energy of − 10.27 kCal/ mol as against binding score energy and hydrogen bond energy of − 41.32 and − 16.31 kCal/mol, respectively, for the novel inhibitor. The molecular dynamics studies also show that all the hydrogen bonds formed by the compound with index number 847b remain stable after the study. 2-(2-(5-amino-2-((((3-aminobenzyl)oxy)carbonyl) amino)-5-oxopentanamido)-4-(2-(tert-butyl)-4-oxo-4-(pentan-3-ylamino) butanamido)-3-hydroxybutyl) benzoic acid differ significantly in activity from its parent chain because of the introduction of primary amine group attached to meta position of first benzene ring and the carboxyl group attached to the ortho position of the second benzene ring. These groups have the ability to increase the overall binding energy by increasing the number of hydrogen bonds interaction present in their complex. This effects makes 2-(2-(5-amino-2-((((3aminobenzyl)oxy)carbonyl) amino)-5-oxopentanamido)-4 -( 2 -( te r t-b u t y l ) -4 -oxo -4 -( p e nt a n -3 -y l a m i n o ) butanamido)-3-hydroxybutyl) benzoic acid a better drug candidate against SARS CoV 2 main protease with the binding energy of -41.32 kCal/mol.
The result could be partly explained by the fact that the inhibitor has a strong hydrogen bond interaction with the amino acids of the binding pocket of the SARS CoV 2 main protease which is evidenced by the high hydrogen bond energy value of − 11.45 kCal/mol making it the highest of all the docked inhibitors. Other noticeable interactions with the receptor include carbon-hydrogen interaction with Thr26 and Thr25, π-sulfur interaction with Met49, amide-pi Stacked interaction with Leu141 and Asn142, π-alkyl interaction with Hie41. Zhijian Xu et al. in their work used both MM/GBSA and SIE methods and they voted for nelfinavir, with the binding free energy of − 24.69 ± 0.52 kCal/mol and − 9.42 ± 0.04 kCal/mol, respectively, to be a potential inhibitor against 2019-nCov Mpro (Xu et al. 2020). The inhibiting capacity of their proposed drug is not comparable to that obtained with the compounds with index numbers 441 and 741, making it a better drug candidate than nelfinavir.
The result for compound with index numbers 441 is shown in Fig. 4, the binding energy is reported in Additional file 1: Table S2 to be − 15.67 kCal/mol and the interaction type result is as shown in Table 4. The docked result shows that the inhibitor has seven hydrogen bond interactions with five amino acids (Thr26, Cys44, Gln189, Thr25, and Met49). The binding energy of inhibitor with Index number 741 is determined to be − 47.88 kCal/mol. This makes it the 2nd most active chemical agent with the ability to inhibit SARS CoV 2 main protease (Additional file 1: Table S2). The docked result owes its binding affinity to the presence of seven H-bond with the amino acids which include Thr26, Glu166, Hie164, and Met165 (Fig. 5). Other interactions such as Π-sulfur-type with Cys145, Met165, π-pi stacked with Hie41 and π-alkyl with Met49 are also observed (Table 5).
Based on hydrogen bond energy ranking, it occupies the third position, and it has hydrogen bond energy of 9.57 kCal/mol. Other important interactions such as π-alkyl, π-sulfur interactions are also reported. This is far better than all the compounds obtained by Motiwale and colleagues (Motiwale et al. 2020). They applied molecular docking approach in conjugation with molecular dynamics (MD) simulations to find out potential inhibitors against Mpro of SARS CoV-2 from previously reported SARS-3CL protease inhibitors. They used a total of 61, previously known inhibitors. According to the molecular docking results, it was found that the binding energy of co-crystallized ligand, JFM (N-(2-phenylethyl)-methanesulfonamide) was found to be − 5.1 kCal/mol while the binding energy of all the 61 inhibitors lies between − 6.2 and − 8.8 kCal/ mol.

SARS CoV 2 main protease (PDB ID: 6XBH) with inhibitor index number 847
As presented in Additional file 1: Table S2, the binding energy of this interaction is reported to be − 39.89 kCal/mol. This is the inhibitor with the highest number of hydrogen bond making it one of the best chemical agents with the ability to inhibit SARS CoV 2 main protease. The docked result owes its binding affinity to the   Other interactions such as π-anion interaction with Glu166 and π-alkyl interaction with Hie41, Met165 also contributed to the high affinity of the inhibitor in the binding site by stabilizing its structure to conform to the surface of the polar amino acids. From the virtual screening results by Khan and colleague, two drug molecules were selected for each drug target protein [Paritaprevir (ΔG = − 9.8 kCal/mol) &Raltegravir (ΔG = − 7.8 kCal/ mol) for 3CLpro and Dolutegravir (ΔG = − 9.4 kCal/mol) and Bictegravir (ΔG = − 8.4 kCal/mol) for 2'-OMTase]. From their extensive computational analysis, they proposed Raltegravir, Paritaprevir, Bictegravir and Dolutegravir as excellent lead candidates for these crucial proteins and they could become potential therapeutic drugs against 2019-nCoV (Khan et al. 2020). This result cannot be compared with our proposed drug that has a binding free energy of − 39.89 kCal/mol.

Conclusions
The molecular docking results shown in the figures confirm that the hydrophobic and hydrogen bonding interactions with these targets have pivotal contributions to the binding structures and binding free energies, even though the van der Waals and π-interactions contributed to the stabilization of the binding structures.
The molecular docking result also shows that, inhibitors with Index numbers 331, 741, 819, 441, 847, and 46 with ICM score of − 48.38 kCal/mol, − 47.88 kCal/mol, − 47.52 kCal/mol, 29.01 kCal/mol, 39.89 kCal/mol, and − 15.67 kCal/mol, respectively, best inhibit SARS CoV 2 main protease of the compounds within our data set. These compounds were further utilized in designing new potent inhibitor compounds by attaching potent fragments to the compounds. Most of the newly designed compounds were reported to be more active than the parent structure. This includes compounds with index number 741a, 847b, and 741d with a binding affinity of − 45.33 kCal/mol, − 41.32 kCal/mol and − 40.12 kCal/ mol, respectively. However, compounds with index numbers 741a and 741b and 46d were not considered to be